/*
  File:    random.cpp
  Purpose: To provide superior random number generators, 
           as compared to the standard rand() function. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		   http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>

#include "random.h"

float (*UniformRandom)(void); // long*);
long random_seed = 0;

void SeedRandom( void )
{ random_seed = -time(NULL); if( !UniformRandom ){ UniformRandom = GreatUniformRandom; } }

// in all of these, amazing > great > good

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define RANDEPS 1.2e-7
#define RNMX (1.0-RANDEPS)

#define IM1 2147483563
#define IM2 2147483399
#define AM_Great (1.0/IM1) 
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NDIV_Great (1+IMM1/NTAB)

// the ran1 in Numerical Recipes in C
// a pretty darned good random number generator
// use this if you expect to use fewer than 100,000,000 calls 
// Must be seeded with a negative number!!

float GoodUniformRandom( void ) // long *idum)
{
 long* idum = &random_seed;

 int j;
 long k;
 static long iy=0;
 static long iv[NTAB];
 float temp;
 if (*idum <= 0 || !iy)
 {
  if( -(*idum) < 1 )
  { *idum=1; }
  else
  { *idum = -(*idum); }
  for( j=NTAB+7 ; j>=0 ; j--)
  {
   k=(*idum)/IQ;
   *idum=IA*(*idum-k*IQ)-IR*k;
   if(*idum < 0)
   { *idum += IM; }
   if(j < NTAB)
   { iv[j] = *idum; }
  }
  iy=iv[0];
 }
 k=(*idum)/IQ;
 *idum=IA*(*idum-k*IQ)-IR*k; 
 if(*idum < 0){ *idum += IM; }
 j=iy/NDIV;
 iy=iv[j]; 
 iv[j] = *idum;
 if( (temp=AM*iy) > RNMX )
 { return RNMX; } 
 else
 { return temp; }
}

// the ran2 in Numerical Recipes in C
// a pretty darned great random number generator
// use this if you expect to use greater than 100,000,000 calls 
// Must be seeded with a negative number!!

float GreatUniformRandom(void)// long *idum)
{
 long* idum = &random_seed;
 
 int j;
 long k;
 static long idum2=123456789;
 static long iy=0;
 static long iv[NTAB];
 float temp;
 if( *idum <= 0 )
 {
  if( -(*idum) < 1 )
  { *idum=1; }
  else 
  { *idum = -(*idum); }
  idum2=(*idum);
  for( j=NTAB+7 ; j>=0 ; j--)
  {
   k=(*idum)/IQ1;
   *idum=IA1*(*idum-k*IQ1)-k*IR1;
   if( *idum < 0 )
   { *idum += IM1; }
   if( j < NTAB )
   { iv[j] = *idum; }
  }
  iy=iv[0];
 }
 k=(*idum)/IQ1;
 *idum=IA1*(*idum-k*IQ1)-k*IR1; 
 if (*idum < 0)
 { *idum += IM1; }
 k=idum2/IQ2;
 idum2=IA2*(idum2-k*IQ2)-k*IR2; 
 if( idum2 < 0 )
 { idum2 += IM2; }
 j=iy/NDIV_Great;
 iy=iv[j]-idum2;
 iv[j] = *idum;
 if( iy < 1 )
 { iy += IMM1; }
 if( (temp=AM_Great*iy) > RNMX )
 { return RNMX; }
 else
 { return temp; }
}

#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)

float KnuthUniformRandom(void)// long *idum)
{
 long* idum = &random_seed;
 
 static int inext,inextp;
 static long ma[56]; 
 static int iff=0;
 long mj,mk;
 int i,ii,k;
 if( *idum < 0 || iff == 0 )
 {
  iff=1;
  mj=labs(MSEED-labs(*idum)); 
  mj %= MBIG;
  ma[55]=mj;
  mk=1;
  for( i=1 ; i<=54 ; i++ )
  {
   ii=(21*i) % 55; 
   ma[ii]=mk; 
   mk=mj-mk;
   if (mk < MZ) mk += MBIG;
   mj=ma[ii];
  }
  for( k=1 ; k<=4 ; k++ )
  {
   for( i=1 ; i<=55 ; i++ )
   {
    ma[i] -= ma[1+(i+30) % 55];
    if( ma[i] < MZ )
    { ma[i] += MBIG; }
   }
  }
  inext=0;
  inextp=31; 
  *idum=1;
 }
 if( ++inext == 56 )
 { inext=1; } 
 if( ++inextp == 56 )
 { inextp=1; }
 mj=ma[inext]-ma[inextp];
 if (mj < MZ)
 { mj += MBIG; }
 ma[inext]=mj; 
 return mj*FAC;
}

// gives normal distribution with mean 0, variance 1;

float NormalRandom(void) // long *idum)
{
 long* idum = &random_seed;
 
 if( !UniformRandom )
 { UniformRandom = GoodUniformRandom; }

// float UniformRandom(long *idum);
 static int iset=0;
 static float gset;
 float fac,rsq,v1,v2;
 if( *idum < 0 )
 { iset=0; } 
 if( iset == 0 )
 {
  do
  {
   v1=2.0*UniformRandom()-1.0; // idum)-1.0;
   v2=2.0*UniformRandom()-1.0; // idum)-1.0;
   rsq=v1*v1+v2*v2; 
  }
  while( rsq >= 1.0 || rsq == 0.0 );
  fac=sqrt(-2.0*log(rsq)/rsq);

  gset=v1*fac;
  iset=1; 
  return v2*fac;
 }
 else
 {
  iset=0; 
  return gset;
 }
}

